home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / hyperg_2F1.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  27.1 KB  |  911 lines
  1. /* specfunc/hyperg_2F1.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_pow_int.h>
  27. #include <gsl/gsl_sf_gamma.h>
  28. #include <gsl/gsl_sf_psi.h>
  29. #include <gsl/gsl_sf_hyperg.h>
  30.  
  31. #include "error.h"
  32.  
  33. #define locEPS (1000.0*GSL_DBL_EPSILON)
  34.  
  35.  
  36. /* Assumes c != negative integer.
  37.  */
  38. static int
  39. hyperg_2F1_series(const double a, const double b, const double c,
  40.                   const double x, 
  41.           gsl_sf_result * result
  42.           )
  43. {
  44.   double sum_pos = 1.0;
  45.   double sum_neg = 0.0;
  46.   double del_pos = 1.0;
  47.   double del_neg = 0.0;
  48.   double del = 1.0;
  49.   double k = 0.0;
  50.   int i = 0;
  51.  
  52.   if(fabs(c) < GSL_DBL_EPSILON) {
  53.     result->val = 0.0; /* FIXME: ?? */
  54.     result->err = 1.0;
  55.     GSL_ERROR ("error", GSL_EDOM);
  56.   }
  57.  
  58.   do {
  59.     if(++i > 30000) {
  60.       result->val  = sum_pos - sum_neg;
  61.       result->err  = del_pos + del_neg;
  62.       result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  63.       result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
  64.       GSL_ERROR ("error", GSL_EMAXITER);
  65.     }
  66.     del *= (a+k)*(b+k) * x / ((c+k) * (k+1.0));  /* Gauss series */
  67.  
  68.     if(del > 0.0) {
  69.       del_pos  =  del;
  70.       sum_pos +=  del;
  71.     }
  72.     else if(del == 0.0) {
  73.       /* Exact termination (a or b was a negative integer).
  74.        */
  75.       del_pos = 0.0;
  76.       del_neg = 0.0;
  77.       break;
  78.     }
  79.     else {
  80.       del_neg  = -del;
  81.       sum_neg -=  del;
  82.     }
  83.  
  84.     k += 1.0;
  85.   } while(fabs((del_pos + del_neg)/(sum_pos-sum_neg)) > GSL_DBL_EPSILON);
  86.  
  87.   result->val  = sum_pos - sum_neg;
  88.   result->err  = del_pos + del_neg;
  89.   result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  90.   result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);
  91.  
  92.   return GSL_SUCCESS;
  93. }
  94.  
  95.  
  96. /* a = aR + i aI, b = aR - i aI */
  97. static
  98. int
  99. hyperg_2F1_conj_series(const double aR, const double aI, const double c,
  100.                        double x,
  101.                gsl_sf_result * result)
  102. {
  103.   if(c == 0.0) {
  104.     result->val = 0.0; /* FIXME: should be Inf */
  105.     result->err = 0.0;
  106.     GSL_ERROR ("error", GSL_EDOM);
  107.   }
  108.   else {
  109.     double sum_pos = 1.0;
  110.     double sum_neg = 0.0;
  111.     double del_pos = 1.0;
  112.     double del_neg = 0.0;
  113.     double del = 1.0;
  114.     double k = 0.0;
  115.     do {
  116.       del *= ((aR+k)*(aR+k) + aI*aI)/((k+1.0)*(c+k)) * x;
  117.  
  118.       if(del >= 0.0) {
  119.         del_pos  =  del;
  120.         sum_pos +=  del;
  121.       }
  122.       else {
  123.         del_neg  = -del;
  124.         sum_neg -=  del;
  125.       }
  126.  
  127.       if(k > 30000) {
  128.         result->val  = sum_pos - sum_neg;
  129.         result->err  = del_pos + del_neg;
  130.         result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  131.         result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
  132.         GSL_ERROR ("error", GSL_EMAXITER);
  133.       }
  134.  
  135.       k += 1.0;
  136.     } while(fabs((del_pos + del_neg)/(sum_pos - sum_neg)) > GSL_DBL_EPSILON);
  137.  
  138.     result->val  = sum_pos - sum_neg;
  139.     result->err  = del_pos + del_neg;
  140.     result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  141.     result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);
  142.  
  143.     return GSL_SUCCESS;
  144.   }
  145. }
  146.  
  147.  
  148. /* Luke's rational approximation. The most accesible
  149.  * discussion is in [Kolbig, CPC 23, 51 (1981)].
  150.  * The convergence is supposedly guaranteed for x < 0.
  151.  * You have to read Luke's books to see this and other
  152.  * results. Unfortunately, the stability is not so
  153.  * clear to me, although it seems very efficient when
  154.  * it works.
  155.  */
  156. static
  157. int
  158. hyperg_2F1_luke(const double a, const double b, const double c,
  159.                 const double xin, 
  160.                 gsl_sf_result * result)
  161. {
  162.   int stat_iter;
  163.   const double RECUR_BIG = 1.0e+50;
  164.   const int nmax = 20000;
  165.   int n = 3;
  166.   const double x  = -xin;
  167.   const double x3 = x*x*x;
  168.   const double t0 = a*b/c;
  169.   const double t1 = (a+1.0)*(b+1.0)/(2.0*c);
  170.   const double t2 = (a+2.0)*(b+2.0)/(2.0*(c+1.0));
  171.   double F = 1.0;
  172.   double prec;
  173.  
  174.   double Bnm3 = 1.0;                                  /* B0 */
  175.   double Bnm2 = 1.0 + t1 * x;                         /* B1 */
  176.   double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x);    /* B2 */
  177.  
  178.   double Anm3 = 1.0;                                                      /* A0 */
  179.   double Anm2 = Bnm2 - t0 * x;                                            /* A1 */
  180.   double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x;   /* A2 */
  181.  
  182.   while(1) {
  183.     double npam1 = n + a - 1;
  184.     double npbm1 = n + b - 1;
  185.     double npcm1 = n + c - 1;
  186.     double npam2 = n + a - 2;
  187.     double npbm2 = n + b - 2;
  188.     double npcm2 = n + c - 2;
  189.     double tnm1  = 2*n - 1;
  190.     double tnm3  = 2*n - 3;
  191.     double tnm5  = 2*n - 5;
  192.     double n2 = n*n;
  193.     double F1 =  (3.0*n2 + (a+b-6)*n + 2 - a*b - 2*(a+b)) / (2*tnm3*npcm1);
  194.     double F2 = -(3.0*n2 - (a+b+6)*n + 2 - a*b)*npam1*npbm1/(4*tnm1*tnm3*npcm2*npcm1);
  195.     double F3 = (npam2*npam1*npbm2*npbm1*(n-a-2)*(n-b-2)) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
  196.     double E  = -npam1*npbm1*(n-c-1) / (2*tnm3*npcm2*npcm1);
  197.  
  198.     double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
  199.     double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
  200.     double r = An/Bn;
  201.  
  202.     prec = fabs((F - r)/F);
  203.     F = r;
  204.  
  205.     if(prec < GSL_DBL_EPSILON || n > nmax) break;
  206.  
  207.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  208.       An   /= RECUR_BIG;
  209.       Bn   /= RECUR_BIG;
  210.       Anm1 /= RECUR_BIG;
  211.       Bnm1 /= RECUR_BIG;
  212.       Anm2 /= RECUR_BIG;
  213.       Bnm2 /= RECUR_BIG;
  214.       Anm3 /= RECUR_BIG;
  215.       Bnm3 /= RECUR_BIG;
  216.     }
  217.     else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
  218.       An   *= RECUR_BIG;
  219.       Bn   *= RECUR_BIG;
  220.       Anm1 *= RECUR_BIG;
  221.       Bnm1 *= RECUR_BIG;
  222.       Anm2 *= RECUR_BIG;
  223.       Bnm2 *= RECUR_BIG;
  224.       Anm3 *= RECUR_BIG;
  225.       Bnm3 *= RECUR_BIG;
  226.     }
  227.  
  228.     n++;
  229.     Bnm3 = Bnm2;
  230.     Bnm2 = Bnm1;
  231.     Bnm1 = Bn;
  232.     Anm3 = Anm2;
  233.     Anm2 = Anm1;
  234.     Anm1 = An;
  235.   }
  236.  
  237.   result->val  = F;
  238.   result->err  = 2.0 * fabs(prec * F);
  239.   result->err += 2.0 * GSL_DBL_EPSILON * (n+1.0) * fabs(F);
  240.  
  241.   /* FIXME: just a hack: there's a lot of shit going on here */
  242.   result->err *= 8.0 * (fabs(a) + fabs(b) + 1.0);
  243.  
  244.   stat_iter = (n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );
  245.  
  246.   return stat_iter;
  247. }
  248.  
  249.  
  250. /* Luke's rational approximation for the
  251.  * case a = aR + i aI, b = aR - i aI.
  252.  */
  253. static
  254. int
  255. hyperg_2F1_conj_luke(const double aR, const double aI, const double c,
  256.                      const double xin, 
  257.                      gsl_sf_result * result)
  258. {
  259.   int stat_iter;
  260.   const double RECUR_BIG = 1.0e+50;
  261.   const int nmax = 10000;
  262.   int n = 3;
  263.   const double x = -xin;
  264.   const double x3 = x*x*x;
  265.   const double atimesb = aR*aR + aI*aI;
  266.   const double apb     = 2.0*aR;
  267.   const double t0 = atimesb/c;
  268.   const double t1 = (atimesb +     apb + 1.0)/(2.0*c);
  269.   const double t2 = (atimesb + 2.0*apb + 4.0)/(2.0*(c+1.0));
  270.   double F = 1.0;
  271.   double prec;
  272.  
  273.   double Bnm3 = 1.0;                                  /* B0 */
  274.   double Bnm2 = 1.0 + t1 * x;                         /* B1 */
  275.   double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x);    /* B2 */
  276.  
  277.   double Anm3 = 1.0;                                                      /* A0 */
  278.   double Anm2 = Bnm2 - t0 * x;                                            /* A1 */
  279.   double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x;   /* A2 */
  280.  
  281.   while(1) {
  282.     double nm1 = n - 1;
  283.     double nm2 = n - 2;
  284.     double npam1_npbm1 = atimesb + nm1*apb + nm1*nm1;
  285.     double npam2_npbm2 = atimesb + nm2*apb + nm2*nm2;
  286.     double npcm1 = nm1 + c;
  287.     double npcm2 = nm2 + c;
  288.     double tnm1  = 2*n - 1;
  289.     double tnm3  = 2*n - 3;
  290.     double tnm5  = 2*n - 5;
  291.     double n2 = n*n;
  292.     double F1 =  (3.0*n2 + (apb-6)*n + 2 - atimesb - 2*apb) / (2*tnm3*npcm1);
  293.     double F2 = -(3.0*n2 - (apb+6)*n + 2 - atimesb)*npam1_npbm1/(4*tnm1*tnm3*npcm2*npcm1);
  294.     double F3 = (npam2_npbm2*npam1_npbm1*(nm2*nm2 - nm2*apb + atimesb)) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
  295.     double E  = -npam1_npbm1*(n-c-1) / (2*tnm3*npcm2*npcm1);
  296.  
  297.     double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
  298.     double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
  299.     double r = An/Bn;
  300.  
  301.     prec = fabs(F - r)/fabs(F);
  302.     F = r;
  303.  
  304.     if(prec < GSL_DBL_EPSILON || n > nmax) break;
  305.  
  306.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  307.       An   /= RECUR_BIG;
  308.       Bn   /= RECUR_BIG;
  309.       Anm1 /= RECUR_BIG;
  310.       Bnm1 /= RECUR_BIG;
  311.       Anm2 /= RECUR_BIG;
  312.       Bnm2 /= RECUR_BIG;
  313.       Anm3 /= RECUR_BIG;
  314.       Bnm3 /= RECUR_BIG;
  315.     }
  316.     else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
  317.       An   *= RECUR_BIG;
  318.       Bn   *= RECUR_BIG;
  319.       Anm1 *= RECUR_BIG;
  320.       Bnm1 *= RECUR_BIG;
  321.       Anm2 *= RECUR_BIG;
  322.       Bnm2 *= RECUR_BIG;
  323.       Anm3 *= RECUR_BIG;
  324.       Bnm3 *= RECUR_BIG;
  325.     }
  326.  
  327.     n++;
  328.     Bnm3 = Bnm2;
  329.     Bnm2 = Bnm1;
  330.     Bnm1 = Bn;
  331.     Anm3 = Anm2;
  332.     Anm2 = Anm1;
  333.     Anm1 = An;
  334.   }
  335.   
  336.   result->val  = F;
  337.   result->err  = 2.0 * fabs(prec * F);
  338.   result->err += 2.0 * GSL_DBL_EPSILON * (n+1.0) * fabs(F);
  339.  
  340.   /* FIXME: see above */
  341.   result->err *= 8.0 * (fabs(aR) + fabs(aI) + 1.0);
  342.  
  343.   stat_iter = (n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );
  344.  
  345.   return stat_iter;
  346. }
  347.  
  348.  
  349. /* Do the reflection described in [Moshier, p. 334].
  350.  * Assumes a,b,c != neg integer.
  351.  */
  352. static
  353. int
  354. hyperg_2F1_reflect(const double a, const double b, const double c,
  355.                    const double x, gsl_sf_result * result)
  356. {
  357.   const double d = c - a - b;
  358.   const int intd  = floor(d+0.5);
  359.   const int d_integer = ( fabs(d - intd) < locEPS );
  360.  
  361.   if(d_integer) {
  362.     const double ln_omx = log(1.0 - x);
  363.     const double ad = fabs(d);
  364.     int stat_F2 = GSL_SUCCESS;
  365.     double sgn_2;
  366.     gsl_sf_result F1;
  367.     gsl_sf_result F2;
  368.     double d1, d2;
  369.     gsl_sf_result lng_c;
  370.     gsl_sf_result lng_ad2;
  371.     gsl_sf_result lng_bd2;
  372.     int stat_c;
  373.     int stat_ad2;
  374.     int stat_bd2;
  375.  
  376.     if(d >= 0.0) {
  377.       d1 = d;
  378.       d2 = 0.0;
  379.     }
  380.     else {
  381.       d1 = 0.0;
  382.       d2 = d;
  383.     }
  384.  
  385.     stat_ad2 = gsl_sf_lngamma_e(a+d2, &lng_ad2);
  386.     stat_bd2 = gsl_sf_lngamma_e(b+d2, &lng_bd2);
  387.     stat_c   = gsl_sf_lngamma_e(c,    &lng_c);
  388.  
  389.     /* Evaluate F1.
  390.      */
  391.     if(ad < GSL_DBL_EPSILON) {
  392.       /* d = 0 */
  393.       F1.val = 0.0;
  394.       F1.err = 0.0;
  395.     }
  396.     else {
  397.       gsl_sf_result lng_ad;
  398.       gsl_sf_result lng_ad1;
  399.       gsl_sf_result lng_bd1;
  400.       int stat_ad  = gsl_sf_lngamma_e(ad,   &lng_ad);
  401.       int stat_ad1 = gsl_sf_lngamma_e(a+d1, &lng_ad1);
  402.       int stat_bd1 = gsl_sf_lngamma_e(b+d1, &lng_bd1);
  403.  
  404.       if(stat_ad1 == GSL_SUCCESS && stat_bd1 == GSL_SUCCESS && stat_ad == GSL_SUCCESS) {
  405.         /* Gamma functions in the denominator are ok.
  406.      * Proceed with evaluation.
  407.      */
  408.     int i;
  409.         double sum1 = 1.0;
  410.         double term = 1.0;
  411.         double ln_pre1_val = lng_ad.val + lng_c.val + d2*ln_omx - lng_ad1.val - lng_bd1.val;
  412.     double ln_pre1_err = lng_ad.err + lng_c.err + lng_ad1.err + lng_bd1.err + GSL_DBL_EPSILON * fabs(ln_pre1_val);
  413.     int stat_e;
  414.  
  415.         /* Do F1 sum.
  416.          */
  417.         for(i=1; i<ad; i++) {
  418.       int j = i-1;
  419.       term *= (a + d2 + j) * (b + d2 + j) / (1.0 + d2 + j) / i * (1.0-x);
  420.           sum1 += term;
  421.         }
  422.         
  423.     stat_e = gsl_sf_exp_mult_err_e(ln_pre1_val, ln_pre1_err,
  424.                                       sum1, GSL_DBL_EPSILON*fabs(sum1),
  425.                       &F1);
  426.     if(stat_e == GSL_EOVRFLW) {
  427.           OVERFLOW_ERROR(result);
  428.         }
  429.       }
  430.       else {
  431.         /* Gamma functions in the denominator were not ok.
  432.      * So the F1 term is zero.
  433.      */
  434.         F1.val = 0.0;
  435.     F1.err = 0.0;
  436.       }
  437.     } /* end F1 evaluation */
  438.  
  439.  
  440.     /* Evaluate F2.
  441.      */
  442.     if(stat_ad2 == GSL_SUCCESS && stat_bd2 == GSL_SUCCESS) {
  443.       /* Gamma functions in the denominator are ok.
  444.        * Proceed with evaluation.
  445.        */
  446.       const int maxiter = 1000;
  447.       int i;
  448.       double psi_1 = -M_EULER;
  449.       gsl_sf_result psi_1pd; 
  450.       gsl_sf_result psi_apd1;
  451.       gsl_sf_result psi_bpd1;
  452.       int stat_1pd  = gsl_sf_psi_e(1.0 + ad, &psi_1pd);
  453.       int stat_apd1 = gsl_sf_psi_e(a + d1,   &psi_apd1);
  454.       int stat_bpd1 = gsl_sf_psi_e(b + d1,   &psi_bpd1);
  455.       int stat_dall = GSL_ERROR_SELECT_3(stat_1pd, stat_apd1, stat_bpd1);
  456.  
  457.       double psi_val = psi_1 + psi_1pd.val - psi_apd1.val - psi_bpd1.val - ln_omx;
  458.       double psi_err = psi_1pd.err + psi_apd1.err + psi_bpd1.err;
  459.       double fact = 1.0;
  460.       double sum2_val = psi_val;
  461.       double sum2_err = psi_err;
  462.       double ln_pre2_val = lng_c.val + d1*ln_omx - lng_ad2.val - lng_bd2.val;
  463.       double ln_pre2_err = lng_c.val + lng_ad2.val + lng_bd2.val + GSL_DBL_EPSILON * fabs(ln_pre2_val);
  464.       int stat_e;
  465.  
  466.       /* Do F2 sum.
  467.        */
  468.       for(i=1; i<maxiter; i++) {
  469.         int j = i-1;
  470.     double term1 = 1.0/(1.0+j)  + 1.0/(1.0+ad+j);
  471.     double term2 = 1.0/(a+d1+j) + 1.0/(b+d1+j);
  472.         psi_val += term1 - term2;
  473.         psi_err += GSL_DBL_EPSILON * (fabs(term1) + fabs(term2));
  474.         fact *= (a+d1+j)*(b+d1+j)/(ad+j)/i * (1.0-x);
  475.         sum2_val += fact * psi_val;
  476.         sum2_err += fabs(fact * psi_err);
  477.       }
  478.  
  479.       if(i == maxiter) stat_F2 = GSL_EMAXITER;
  480.  
  481.       if(sum2_val == 0.0) {
  482.         F2.val = 0.0;
  483.     F2.err = 0.0;
  484.       }
  485.       else {
  486.         stat_e = gsl_sf_exp_mult_err_e(ln_pre2_val, ln_pre2_err,
  487.                                       sum2_val, sum2_err,
  488.                       &F2);
  489.         if(stat_e == GSL_EOVRFLW) {
  490.           result->val = 0.0;
  491.       result->err = 0.0;
  492.       GSL_ERROR ("error", GSL_EOVRFLW);
  493.         }
  494.       }
  495.       stat_F2 = GSL_ERROR_SELECT_2(stat_F2, stat_dall);
  496.     }
  497.     else {
  498.       /* Gamma functions in the denominator not ok.
  499.        * So the F2 term is zero.
  500.        */
  501.       F2.val = 0.0;
  502.       F2.err = 0.0;
  503.     } /* end F2 evaluation */
  504.  
  505.     sgn_2 = ( GSL_IS_ODD(intd) ? -1.0 : 1.0 );
  506.     result->val  = F1.val + sgn_2 * F2.val;
  507.     result->err  = F1.err + F2. err;
  508.     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(F1.val) + fabs(F2.val));
  509.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  510.     return stat_F2;
  511.   }
  512.   else {
  513.     /* d not an integer */
  514.  
  515.     gsl_sf_result pre1, pre2;
  516.     double sgn1, sgn2;
  517.     gsl_sf_result F1, F2;
  518.     int status_F1, status_F2;
  519.  
  520.     /* These gamma functions appear in the denominator, so we
  521.      * catch their harmless domain errors and set the terms to zero.
  522.      */
  523.     gsl_sf_result ln_g1ca,  ln_g1cb,  ln_g2a,  ln_g2b;
  524.     double sgn_g1ca, sgn_g1cb, sgn_g2a, sgn_g2b;
  525.     int stat_1ca = gsl_sf_lngamma_sgn_e(c-a, &ln_g1ca, &sgn_g1ca);
  526.     int stat_1cb = gsl_sf_lngamma_sgn_e(c-b, &ln_g1cb, &sgn_g1cb);
  527.     int stat_2a  = gsl_sf_lngamma_sgn_e(a, &ln_g2a, &sgn_g2a);
  528.     int stat_2b  = gsl_sf_lngamma_sgn_e(b, &ln_g2b, &sgn_g2b);
  529.     int ok1 = (stat_1ca == GSL_SUCCESS && stat_1cb == GSL_SUCCESS);
  530.     int ok2 = (stat_2a  == GSL_SUCCESS && stat_2b  == GSL_SUCCESS);
  531.     
  532.     gsl_sf_result ln_gc,  ln_gd,  ln_gmd;
  533.     double sgn_gc, sgn_gd, sgn_gmd;
  534.     gsl_sf_lngamma_sgn_e( c, &ln_gc,  &sgn_gc);
  535.     gsl_sf_lngamma_sgn_e( d, &ln_gd,  &sgn_gd);
  536.     gsl_sf_lngamma_sgn_e(-d, &ln_gmd, &sgn_gmd);
  537.     
  538.     sgn1 = sgn_gc * sgn_gd  * sgn_g1ca * sgn_g1cb;
  539.     sgn2 = sgn_gc * sgn_gmd * sgn_g2a  * sgn_g2b;
  540.  
  541.     if(ok1 && ok2) {
  542.       double ln_pre1_val = ln_gc.val + ln_gd.val  - ln_g1ca.val - ln_g1cb.val;
  543.       double ln_pre2_val = ln_gc.val + ln_gmd.val - ln_g2a.val  - ln_g2b.val + d*log(1.0-x);
  544.       double ln_pre1_err = ln_gc.err + ln_gd.err + ln_g1ca.err + ln_g1cb.err;
  545.       double ln_pre2_err = ln_gc.err + ln_gmd.err + ln_g2a.err  + ln_g2b.err;
  546.       if(ln_pre1_val < GSL_LOG_DBL_MAX && ln_pre2_val < GSL_LOG_DBL_MAX) {
  547.         gsl_sf_exp_err_e(ln_pre1_val, ln_pre1_err, &pre1);
  548.     gsl_sf_exp_err_e(ln_pre2_val, ln_pre2_err, &pre2);
  549.         pre1.val *= sgn1;
  550.         pre2.val *= sgn2;
  551.       }
  552.       else {
  553.         OVERFLOW_ERROR(result);
  554.       }
  555.     }
  556.     else if(ok1 && !ok2) {
  557.       double ln_pre1_val = ln_gc.val + ln_gd.val - ln_g1ca.val - ln_g1cb.val;
  558.       double ln_pre1_err = ln_gc.err + ln_gd.err + ln_g1ca.err + ln_g1cb.err;
  559.       if(ln_pre1_val < GSL_LOG_DBL_MAX) {
  560.         gsl_sf_exp_err_e(ln_pre1_val, ln_pre1_err, &pre1);
  561.         pre1.val *= sgn1;
  562.         pre2.val = 0.0;
  563.     pre2.err = 0.0;
  564.       }
  565.       else {
  566.         OVERFLOW_ERROR(result);
  567.       }
  568.     }
  569.     else if(!ok1 && ok2) {
  570.       double ln_pre2_val = ln_gc.val + ln_gmd.val - ln_g2a.val - ln_g2b.val + d*log(1.0-x);
  571.       double ln_pre2_err = ln_gc.err + ln_gmd.err + ln_g2a.err + ln_g2b.err;
  572.       if(ln_pre2_val < GSL_LOG_DBL_MAX) {
  573.         pre1.val = 0.0;
  574.         pre1.err = 0.0;
  575.         gsl_sf_exp_err_e(ln_pre2_val, ln_pre2_err, &pre2);
  576.         pre2.val *= sgn2;
  577.       }
  578.       else {
  579.         OVERFLOW_ERROR(result);
  580.       }
  581.     }
  582.     else {
  583.       pre1.val = 0.0;
  584.       pre2.val = 0.0;
  585.       UNDERFLOW_ERROR(result);
  586.     }
  587.  
  588.     status_F1 = hyperg_2F1_series(  a,   b, 1.0-d, 1.0-x, &F1);
  589.     status_F2 = hyperg_2F1_series(c-a, c-b, 1.0+d, 1.0-x, &F2);
  590.  
  591.     result->val  = pre1.val*F1.val + pre2.val*F2.val;
  592.     result->err  = fabs(pre1.val*F1.err) + fabs(pre2.val*F2.err);
  593.     result->err += fabs(pre1.err*F1.val) + fabs(pre2.err*F2.val);
  594.     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(pre1.val*F1.val) + fabs(pre2.val*F2.val));
  595.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  596.  
  597.     return GSL_SUCCESS;
  598.   }
  599. }
  600.  
  601.  
  602. static int pow_omx(const double x, const double p, gsl_sf_result * result)
  603. {
  604.   double ln_omx;
  605.   double ln_result;
  606.   if(fabs(x) < GSL_ROOT5_DBL_EPSILON) {
  607.     ln_omx = -x*(1.0 + x*(1.0/2.0 + x*(1.0/3.0 + x/4.0 + x*x/5.0)));
  608.   }
  609.   else {
  610.     ln_omx = log(1.0-x);
  611.   }
  612.   ln_result = p * ln_omx;
  613.   return gsl_sf_exp_err_e(ln_result, GSL_DBL_EPSILON * fabs(ln_result), result);
  614. }
  615.  
  616.  
  617. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  618.  
  619. int
  620. gsl_sf_hyperg_2F1_e(double a, double b, const double c,
  621.                        const double x,
  622.                        gsl_sf_result * result)
  623. {
  624.   const double d = c - a - b;
  625.   const double rinta = floor(a + 0.5);
  626.   const double rintb = floor(b + 0.5);
  627.   const double rintc = floor(c + 0.5);
  628.   const int a_neg_integer = ( a < 0.0  &&  fabs(a - rinta) < locEPS );
  629.   const int b_neg_integer = ( b < 0.0  &&  fabs(b - rintb) < locEPS );
  630.   const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );
  631.  
  632.   result->val = 0.0;
  633.   result->err = 0.0;
  634.  
  635.   if(x < -1.0 || 1.0 <= x) {
  636.     DOMAIN_ERROR(result);
  637.   }
  638.  
  639.   if(c_neg_integer) {
  640.     if(! (a_neg_integer && a > c + 0.1)) DOMAIN_ERROR(result);
  641.     if(! (b_neg_integer && b > c + 0.1)) DOMAIN_ERROR(result);
  642.   }
  643.  
  644.   if(fabs(c-b) < locEPS || fabs(c-a) < locEPS) {
  645.     return pow_omx(x, d, result);  /* (1-x)^(c-a-b) */
  646.   }
  647.  
  648.   if(a >= 0.0 && b >= 0.0 && c >=0.0 && x >= 0.0) {
  649.     /* Series has all positive definite terms.
  650.      */
  651.     return hyperg_2F1_series(a, b, c, x, result);
  652.   }
  653.  
  654.   if(fabs(a) < 10.0 && fabs(b) < 10.0) {
  655.     /* a and b are not too large, so we attempt
  656.      * variations on the series summation.
  657.      */
  658.     if(a_neg_integer) {
  659.       return hyperg_2F1_series(rinta, b, c, x, result);
  660.     }
  661.     if(b_neg_integer) {
  662.       return hyperg_2F1_series(a, rintb, c, x, result);
  663.     }
  664.  
  665.     if(x < -0.25) {
  666.       return hyperg_2F1_luke(a, b, c, x, result);
  667.     }
  668.     else if(x < 0.5) {
  669.       return hyperg_2F1_series(a, b, c, x, result);
  670.     }
  671.     else {
  672.       if(fabs(c) > 10.0) {
  673.         return hyperg_2F1_series(a, b, c, x, result);
  674.       }
  675.       else {
  676.         return hyperg_2F1_reflect(a, b, c, x, result);
  677.       }
  678.     }
  679.   }
  680.   else {
  681.     /* Either a or b or both large.
  682.      * Introduce some new variables ap,bp so that bp is
  683.      * the larger in magnitude.
  684.      */
  685.     double ap, bp; 
  686.     if(fabs(a) > fabs(b)) {
  687.       bp = a;
  688.       ap = b;
  689.     }
  690.     else {
  691.       bp = b;
  692.       ap = a;
  693.     }
  694.  
  695.     if(x < 0.0) {
  696.       /* What the hell, maybe Luke will converge.
  697.        */
  698.       return hyperg_2F1_luke(a, b, c, x, result);
  699.     }
  700.  
  701.     if(GSL_MAX_DBL(fabs(a),1.0)*fabs(bp)*fabs(x) < 2.0*fabs(c)) {
  702.       /* If c is large enough or x is small enough,
  703.        * we can attempt the series anyway.
  704.        */
  705.       return hyperg_2F1_series(a, b, c, x, result);
  706.     }
  707.  
  708.     if(fabs(bp*bp*x*x) < 0.001*fabs(bp) && fabs(a) < 10.0) {
  709.       /* The famous but nearly worthless "large b" asymptotic.
  710.        */
  711.       int stat = gsl_sf_hyperg_1F1_e(a, c, bp*x, result);
  712.       result->err = 0.001 * fabs(result->val);
  713.       return stat;
  714.     }
  715.  
  716.     /* We give up. */
  717.     result->val = 0.0;
  718.     result->err = 0.0;
  719.     GSL_ERROR ("error", GSL_EUNIMPL);
  720.   }
  721. }
  722.  
  723.  
  724. int
  725. gsl_sf_hyperg_2F1_conj_e(const double aR, const double aI, const double c,
  726.                             const double x,
  727.                             gsl_sf_result * result)
  728. {
  729.   const double ax = fabs(x);
  730.   const double rintc = floor(c + 0.5);
  731.   const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );
  732.  
  733.   result->val = 0.0;
  734.   result->err = 0.0;
  735.  
  736.   if(ax >= 1.0 || c_neg_integer || c == 0.0) {
  737.     DOMAIN_ERROR(result);
  738.   }
  739.  
  740.   if(   (ax < 0.25 && fabs(aR) < 20.0 && fabs(aI) < 20.0)
  741.      || (c > 0.0 && x > 0.0)
  742.     ) {
  743.     return hyperg_2F1_conj_series(aR, aI, c, x, result);
  744.   }
  745.   else if(fabs(aR) < 10.0 && fabs(aI) < 10.0) {
  746.     if(x < -0.25) {
  747.       return hyperg_2F1_conj_luke(aR, aI, c, x, result);
  748.     }
  749.     else {
  750.       return hyperg_2F1_conj_series(aR, aI, c, x, result);
  751.     }
  752.   }
  753.   else {
  754.     if(x < 0.0) {
  755.       /* What the hell, maybe Luke will converge.
  756.        */
  757.       return hyperg_2F1_conj_luke(aR, aI, c, x, result); 
  758.     }
  759.  
  760.     /* Give up. */
  761.     result->val = 0.0;
  762.     result->err = 0.0;
  763.     GSL_ERROR ("error", GSL_EUNIMPL);
  764.   }
  765. }
  766.  
  767.  
  768. int
  769. gsl_sf_hyperg_2F1_renorm_e(const double a, const double b, const double c,
  770.                               const double x,
  771.                   gsl_sf_result * result
  772.                   )
  773. {
  774.   const double rinta = floor(a + 0.5);
  775.   const double rintb = floor(b + 0.5);
  776.   const double rintc = floor(c + 0.5);
  777.   const int a_neg_integer = ( a < 0.0  &&  fabs(a - rinta) < locEPS );
  778.   const int b_neg_integer = ( b < 0.0  &&  fabs(b - rintb) < locEPS );
  779.   const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );
  780.   
  781.   if(c_neg_integer) {
  782.     if((a_neg_integer && a > c+0.1) || (b_neg_integer && b > c+0.1)) {
  783.       /* 2F1 terminates early */
  784.       result->val = 0.0;
  785.       result->err = 0.0;
  786.       return GSL_SUCCESS;
  787.     }
  788.     else {
  789.       /* 2F1 does not terminate early enough, so something survives */
  790.       /* [Abramowitz+Stegun, 15.1.2] */
  791.       gsl_sf_result g1, g2, g3, g4, g5;
  792.       double s1, s2, s3, s4, s5;
  793.       int stat = 0;
  794.       stat += gsl_sf_lngamma_sgn_e(a-c+1, &g1, &s1);
  795.       stat += gsl_sf_lngamma_sgn_e(b-c+1, &g2, &s2);
  796.       stat += gsl_sf_lngamma_sgn_e(a, &g3, &s3);
  797.       stat += gsl_sf_lngamma_sgn_e(b, &g4, &s4);
  798.       stat += gsl_sf_lngamma_sgn_e(-c+2, &g5, &s5);
  799.       if(stat != 0) {
  800.         DOMAIN_ERROR(result);
  801.       }
  802.       else {
  803.         gsl_sf_result F;
  804.         int stat_F = gsl_sf_hyperg_2F1_e(a-c+1, b-c+1, -c+2, x, &F);
  805.         double ln_pre_val = g1.val + g2.val - g3.val - g4.val - g5.val;
  806.     double ln_pre_err = g1.err + g2.err + g3.err + g4.err + g5.err;
  807.     double sg  = s1 * s2 * s3 * s4 * s5;
  808.     int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
  809.                                           sg * F.val, F.err,
  810.                           result);
  811.     return GSL_ERROR_SELECT_2(stat_e, stat_F);
  812.       }
  813.     }
  814.   }
  815.   else {
  816.     /* generic c */
  817.     gsl_sf_result F;
  818.     gsl_sf_result lng;
  819.     double sgn;
  820.     int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);
  821.     int stat_F = gsl_sf_hyperg_2F1_e(a, b, c, x, &F);
  822.     int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,
  823.                                           sgn*F.val, F.err,
  824.                       result);
  825.     return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);
  826.   }
  827. }
  828.  
  829.  
  830. int
  831. gsl_sf_hyperg_2F1_conj_renorm_e(const double aR, const double aI, const double c,
  832.                                    const double x,
  833.                        gsl_sf_result * result
  834.                        )
  835. {
  836.   const double rintc = floor(c  + 0.5);
  837.   const double rinta = floor(aR + 0.5);
  838.   const int a_neg_integer = ( aR < 0.0 && fabs(aR-rinta) < locEPS && aI == 0.0);
  839.   const int c_neg_integer = (  c < 0.0 && fabs(c - rintc) < locEPS );
  840.  
  841.   if(c_neg_integer) {
  842.     if(a_neg_integer && aR > c+0.1) {
  843.       /* 2F1 terminates early */
  844.       result->val = 0.0;
  845.       result->err = 0.0;
  846.       return GSL_SUCCESS;
  847.     }
  848.     else {
  849.       /* 2F1 does not terminate early enough, so something survives */
  850.       /* [Abramowitz+Stegun, 15.1.2] */
  851.       gsl_sf_result g1, g2;
  852.       gsl_sf_result g3;
  853.       gsl_sf_result a1, a2;
  854.       int stat = 0;
  855.       stat += gsl_sf_lngamma_complex_e(aR-c+1, aI, &g1, &a1);
  856.       stat += gsl_sf_lngamma_complex_e(aR, aI, &g2, &a2);
  857.       stat += gsl_sf_lngamma_e(-c+2.0, &g3);
  858.       if(stat != 0) {
  859.         DOMAIN_ERROR(result);
  860.       }
  861.       else {
  862.         gsl_sf_result F;
  863.         int stat_F = gsl_sf_hyperg_2F1_conj_e(aR-c+1, aI, -c+2, x, &F);
  864.         double ln_pre_val = 2.0*(g1.val - g2.val) - g3.val;
  865.     double ln_pre_err = 2.0 * (g1.err + g2.err) + g3.err;
  866.     int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
  867.                                               F.val, F.err,
  868.                                               result);
  869.     return GSL_ERROR_SELECT_2(stat_e, stat_F);
  870.       }
  871.     }
  872.   }
  873.   else {
  874.     /* generic c */
  875.     gsl_sf_result F;
  876.     gsl_sf_result lng;
  877.     double sgn;
  878.     int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);
  879.     int stat_F = gsl_sf_hyperg_2F1_conj_e(aR, aI, c, x, &F);
  880.     int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,
  881.                                           sgn*F.val, F.err,
  882.                                           result);
  883.     return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);
  884.   }
  885. }
  886.  
  887.  
  888. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  889.  
  890. #include "eval.h"
  891.  
  892. double gsl_sf_hyperg_2F1(double a, double b, double c, double x)
  893. {
  894.   EVAL_RESULT(gsl_sf_hyperg_2F1_e(a, b, c, x, &result));
  895. }
  896.  
  897. double gsl_sf_hyperg_2F1_conj(double aR, double aI, double c, double x)
  898. {
  899.   EVAL_RESULT(gsl_sf_hyperg_2F1_conj_e(aR, aI, c, x, &result));
  900. }
  901.  
  902. double gsl_sf_hyperg_2F1_renorm(double a, double b, double c, double x)
  903. {
  904.   EVAL_RESULT(gsl_sf_hyperg_2F1_renorm_e(a, b, c, x, &result));
  905. }
  906.  
  907. double gsl_sf_hyperg_2F1_conj_renorm(double aR, double aI, double c, double x)
  908. {
  909.   EVAL_RESULT(gsl_sf_hyperg_2F1_conj_renorm_e(aR, aI, c, x, &result));
  910. }
  911.